function [is_near ds_near]= SamplesAdjacency(x_n,x_s,rho)
% function [is_near ds_near]= SamplesAdjacency(x_n,x_s,range)
% This funtion computes to each sample point its neihbors in x_n

range_n = 2*rho;
n_n = size(x_n,1);
n_s = size(x_s,1);
dim = size(x_n,2);

if n_s==1
    x_s=[x_s;x_s];
    n_s=2;
end

if dim>1
  atria  = nn_prepare(x_s);
end

if dim>n_n
  error('DIM is bigger than number of nodes, maybe you needs transpose your data')
end
is_near  = {zeros(n_s,1)};
ds_near  = {zeros(n_s,1)};
for k=1:n_s
  is_near{k} = [];
  ds_near{k} = [];
end

for i=1:n_n
  if dim>1
    [count, nd_near] = range_search(x_s, atria, x_n(i,:), range_n(i));
  else
    [count, nd_near] = range_search1D(x_s, x_n(i), range_n(i));
  end
  if count > 0
    %     n_near = sort(n_near{1,1});
    n_near = nd_near{1,1};
    d_near = nd_near{1,2};
    for j=1:count
      s   = n_near(j);
      d   = d_near(j);
      is_near{s} =  [is_near{s} i];
      ds_near{s} =  [ds_near{s} d];
    end
  end

end


for k=1:n_s
  % the adjacency list is rearanged such that the first indexes are related
  % with the closest nodes to each sample point.
  [dist, ids] = sort(ds_near{k},'ascend');
  is_near{k}  = is_near{k}(ids);
  ds_near{k}  = dist;
end

function [count near] = range_search1D(x_s,x_a,range_a)
% This funtion computes to each x_a node point its neihbors in x_s

% Naive construction of the neighbor list
sPts = length(x_s);
all_ = (1:sPts);
% Naive find in range
dist = (x_s-x_a).^2;
dist = sqrt(dist);

near_= all_(dist<range_a);
count= length(near_);
dist_= dist(near_);
near = {near_,dist_};
